Data Preparation
#For forecasting, we chose the latest data
trend_death_hb <- trend_hb_daily %>%
filter (hb_name == "Scotland") %>%
filter(date >="2021-06-01") %>%
filter(!(is.na(daily_deaths))) %>%
select(date, daily_deaths)
# Convert it into a time series
daily_death_hb_zoo <- zoo(trend_death_hb$daily_deaths,
order.by=as.Date(trend_death_hb$date, format='%m/%d/%Y'))
# Convert it into a time series
daily_death_hb_timeseries <- timeSeries::as.timeSeries(daily_death_hb_zoo)
Step 1 : Visualize the time series
original_series_death<-autoplot(daily_death_hb_timeseries, ts.colour = '#5ab4ac')+
xlab("Month") +
ylab("Patient died")+
#scale_x_date(breaks = "1 month", date_labels = "%b - %y" )+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
ggtitle("Trend on Deaths") +
color_theme()
ggplotly(original_series_death)
Step 2 : Identification of model : (Finding d:)
Identify whether the time series is stationary / non stationary we can use ADF Augmented Dickey-Fuller test
adf_test_death <- adf.test(daily_death_hb_timeseries)
adf_test_death
Augmented Dickey-Fuller Test
data: daily_death_hb_timeseries
Dickey-Fuller = -1.4701, Lag order = 4, p-value = 0.7967
alternative hypothesis: stationary
The time series is not stationary since we have a high p-value (p-value must be < 0.05). So we apply difference
first_diff_death<- diff(daily_death_hb_timeseries)
adf_test1_death <- adf.test(na.omit(first_diff_death))
Warning in adf.test(na.omit(first_diff_death)) :
p-value smaller than printed p-value
adf_test1_death
Augmented Dickey-Fuller Test
data: na.omit(first_diff_death)
Dickey-Fuller = -5.5546, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary
Create a dataframe to compare
adf_data_death <- data.frame(Data = c("Original", "First-Ordered"),
Dickey_Fuller = c(adf_test_death$statistic, adf_test1_death$statistic),
p_value = c(adf_test_death$p.value,adf_test1_death$p.value))
adf_data_death
Initially the p-value is high which indicates that the Time Series is not stationary. So we apply difference 1 time. After the first difference, the p-value < significance level (0.05) So we can conclude that the difference data are stationary. So difference (d = 1)
Other method:
ndiffs(daily_death_hb_timeseries)
Let’s plot the First Order Difference Series
Order of first difference
first_diff_death<- diff(daily_death_hb_timeseries)
p<- autoplot(first_diff_death, ts.colour = '#5ab4ac') +
xlab("Month") +
ylab("DEATH")+
# scale_x_date(breaks = "1 month", date_labels = "%b - %y" )+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
ggtitle("First-Order Difference Series") +
color_theme()
ggplotly(p)
Step 3 Estimate the parameters (Finding p and q)
For our model ARIMA (p,d,q), we found d = 1, the next step is to get the values of p and q, the order of AR and MA part. Plot ACF and PACF charts to identify q and p respectively.
par(mfrow=c(2,2))
acf_death <- acf1(first_diff_death, col=2:7, lwd=4)
pacf_death <- acf1(first_diff_death, pacf = TRUE, col=2:7, lwd=4)
The ACF and PACF plots tapers to zeo at one point. So it follows an ARMA model
The PACF non-zero value at lag 3 So the data may follow an AR(3) model
The ACF non - zero value at lag 2 So it can be MA (2)
So we propose three ARMA models for the differenced data: ARMA(p,q) ARMA(3,2), ARMA(3,0) and ARMA(0,2).
That is, for the original time series, we propose three ARIMA models,ARIMA(p,d,q) ARIMA(3,1,2), ARIMA(3,1,0) and ARMA(0,1,2).
Step 4 Build the ARIMA model
Manual ARIMA
arima_fit_dth_1 = Arima(daily_death_hb_timeseries, order = c(3,1,2))
arima_fit_dth_2 = Arima(daily_death_hb_timeseries, order = c(3,1,0))
arima_fit_dth_3 = Arima(daily_death_hb_timeseries, order = c(0,1,2))
summary(arima_fit_dth_1)
summary(arima_fit_dth_2)
summary(arima_fit_dth_3)
Another way of checking AIC
texreg::screenreg(list(arima_fit_dth_1, arima_fit_dth_2, arima_fit_dth_3),
custom.model.names =c("ARIMA(3,1,2)","ARIMA(3,1,0)","ARIMA(0,1,2)"),
center = TRUE,
table = FALSE)
Based on this, ARIMA model (3,1,2) seems best. We can verify the same using automated method
Automated ARIMA
auto_arima_fit_death <- auto.arima(lag(daily_death_hb_timeseries),
seasonal=FALSE,
stepwise=FALSE,
approximation=FALSE,
trace = TRUE
)
ARIMA(0,1,0) : 661.2495
ARIMA(0,1,0) with drift : 662.9805
ARIMA(0,1,1) : 604.9241
ARIMA(0,1,1) with drift : 603.4592
ARIMA(0,1,2) : 600.5927
ARIMA(0,1,2) with drift : 600.1606
ARIMA(0,1,3) : 602.1154
ARIMA(0,1,3) with drift : 601.3862
ARIMA(0,1,4) : 597.6822
ARIMA(0,1,4) with drift : 597.7631
ARIMA(0,1,5) : 599.8704
ARIMA(0,1,5) with drift : 599.9091
ARIMA(1,1,0) : 605.4993
ARIMA(1,1,0) with drift : 606.5658
ARIMA(1,1,1) : 599.2604
ARIMA(1,1,1) with drift : 598.7763
ARIMA(1,1,2) : 601.2259
ARIMA(1,1,2) with drift : 600.6498
ARIMA(1,1,3) : 603.1702
ARIMA(1,1,3) with drift : 602.3916
ARIMA(1,1,4) : 599.8811
ARIMA(1,1,4) with drift : 599.963
ARIMA(2,1,0) : 601.7416
ARIMA(2,1,0) with drift : 602.2924
ARIMA(2,1,1) : 601.1992
ARIMA(2,1,1) with drift : 600.5597
ARIMA(2,1,2) : 603.3692
ARIMA(2,1,2) with drift : 602.7661
ARIMA(2,1,3) : 592.0755
ARIMA(2,1,3) with drift : Inf
ARIMA(3,1,0) : 603.1214
ARIMA(3,1,0) with drift : 603.4379
ARIMA(3,1,1) : 603.3196
ARIMA(3,1,1) with drift : 602.7205
ARIMA(3,1,2) : 599.3322
ARIMA(3,1,2) with drift : 600.2376
ARIMA(4,1,0) : 596.0037
ARIMA(4,1,0) with drift : 594.9402
ARIMA(4,1,1) : 594.7462
ARIMA(4,1,1) with drift : 594.1175
ARIMA(5,1,0) : 593.284
ARIMA(5,1,0) with drift : 593.3211
Best model: ARIMA(2,1,3)
auto_arima_fit_death
Series: lag(daily_death_hb_timeseries)
ARIMA(2,1,3)
Coefficients:
ar1 ar2 ma1 ma2 ma3
-1.6080 -0.9445 0.8866 0.0474 -0.4563
s.e. 0.0583 0.0518 0.0987 0.1133 0.0984
sigma^2 estimated as 7.525: log likelihood=-289.67
AIC=591.33 AICc=592.08 BIC=608.06
However Automated ARIMA also confirms that the ARIMA(2,1,3) seems good based on AIC
coef_dth<-lmtest::coeftest(auto_arima_fit_death)
coef_dth
Model Selection Criteria :
ARIMA models with minimum AIC, RMSE and MAPE criteria were chosen as the best models. Based on Akaike Information Criterion (AIC) above, an ARIMA(2, 1, 3) model seems best.
Step 5 Check for Diagnostics
Let’s plot the diagnostics with the results to make sure the normality and correlation assumptions for the model hold. If the residuals look like white noise, proceed with forecast and prediction, otherwise repeat the model building.
res_dth <-checkresiduals(auto_arima_fit_death, theme = color_theme())
res_dth
The ACF plot of the residuals from the ARIMA(2,1,3) model shows that all auto correlations are within the threshold limits except 1, But portmanteau test (Ljung-Box test) shows a higher p-value which means the values are indpendent. i.e We fail to reject the null hypothesis
Fitting the ARIMA model with the existing data
Let’s plot the actuals against the fitted values
Convert model and time series to dataframe for plotting
daily_death_hb_timeseries_data <- fortify(daily_death_hb_timeseries) %>%
clean_names() %>%
remove_rownames %>%
rename (date = index,
death = data)%>%
mutate(index = seq(1:nrow(daily_death_hb_timeseries)))
arima_fit_dth_resid <- ts(daily_death_hb_timeseries[1:nrow(daily_death_hb_timeseries)]) - resid(auto_arima_fit_death)
arima_fit_dth_data <- fortify(arima_fit_dth_resid) %>%
clean_names() %>%
mutate(data = round(data,2))%>%
mutate (data = ifelse(data < 0,0,data))
fit_existing_dth_data <- daily_death_hb_timeseries_data %>%
inner_join(arima_fit_dth_data, by = c("index"))
Plotting the series along with the fitted values
fit_existing_dth_plot <- fit_existing_dth_data %>%
mutate (data = ifelse(data < 0,0,data)) %>%
mutate(date = as.Date(date)) %>%
ggplot()+
aes(x=date, y = death)+
geom_line(color ="#5ab4ac")+
geom_line(aes(x= date, y = data), colour = "red" )+
xlab("Month") +
ylab("Deaths reported")+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
scale_x_date(breaks = "1 month", date_labels = "%b - %y" )+
ggtitle("Fitting the ARIMA model with existing data") +
color_theme()
ggplotly(fit_existing_dth_plot)
Step 6 Forecast using the model
Data Preparation :
forecast_dth_model <- forecast(auto_arima_fit_death,level = c(80, 95), h = 30)
#Convert the model to dataframe for plotting
# Negative values of the CI interval are considered as 0
forecast_dth_model_data <- fortify(forecast_dth_model) %>%
clean_names() %>%
mutate(data = round(data,2),
fitted= round(fitted,2)) %>%
mutate (lo_80 = ifelse(lo_80 < 0,0,lo_80),
lo_95 = ifelse(lo_95 < 0,0,lo_95)
)
forecast_start_date <- as.Date(max(daily_death_hb_timeseries_data$date)+1)
forecast_end_date <- as.Date(forecast_start_date+29)
forecast_dth_data <- forecast_dth_model_data %>%
filter(!(is.na(point_forecast))) %>%
mutate(date = seq(forecast_start_date,forecast_end_date, by =1)) %>%
select(-data,-fitted, -index)
fitted_dth_data <- forecast_dth_model_data %>%
filter(!(is.na(data))) %>%
inner_join(daily_death_hb_timeseries_data, by = c("index")) %>%
mutate(date = as.Date(date)) %>%
select(date, data, fitted)
Plotting the Vaccination series plus the forecast and 80 - 95% prediction intervals
annotation <- data.frame(
x = c(as.Date("13-06-2021","%d-%m-%Y"),as.Date("11-10-2021","%d-%m-%Y")),
y = c(30,40),
label = c("PAST", "FUTURE")
)
#Time series plots for the next 60 days according to best ARIMA models with 80%–95% CI.
forecast_data_dth_plot <- fitted_dth_data %>%
ggplot()+
geom_line(aes(x= date, y = data), color = "#5ab4ac")+
#geom_line(aes(x= date, y = fitted), colour = "red" )+
geom_line(aes(x= date, y =point_forecast), color ="blue", data = forecast_dth_data )+
geom_ribbon(aes(x = date, y = point_forecast, ymin = lo_80, ymax = hi_80),
data = forecast_dth_data, alpha = 0.3, fill = "green")+
geom_ribbon(aes(x = date, y = point_forecast, ymin = lo_95, ymax = hi_95),
data = forecast_dth_data, alpha = 0.1)+
geom_vline(aes(xintercept=as.numeric(max(date))),color="#f1a340", linetype="dashed",data = fitted_dth_data)+
ggtitle("Projection of new Deaths") +
xlab("Month") +
ylab("Death reported")+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
color_theme()+
scale_x_date(breaks = "1 month", date_labels = "%b - %y" )+
geom_text(data=annotation,
aes( x=x, y=y, label=label),
color="blue",
size=4 )
ggplotly(forecast_data_dth_plot )
LS0tDQp0aXRsZTogIkZvcmVjYXN0IG9uIERlYXRocyAoQVJJTUEgTW9kZWxsaW5nKSINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KZWRpdG9yX29wdGlvbnM6IA0KICBtYXJrZG93bjogDQogICAgd3JhcDogc2VudGVuY2UNCi0tLQ0KDQoqKkRhdGEgUHJlcGFyYXRpb24qKg0KDQpgYGB7cn0NCiNGb3IgZm9yZWNhc3RpbmcsIHdlIGNob3NlIHRoZSBsYXRlc3QgZGF0YQ0KdHJlbmRfZGVhdGhfaGIgPC0gdHJlbmRfaGJfZGFpbHkgJT4lIA0KICBmaWx0ZXIgKGhiX25hbWUgPT0gIlNjb3RsYW5kIikgJT4lIA0KICBmaWx0ZXIoZGF0ZSA+PSIyMDIxLTA2LTAxIikgJT4lIA0KICBmaWx0ZXIoIShpcy5uYShkYWlseV9kZWF0aHMpKSkgJT4lIA0KICBzZWxlY3QoZGF0ZSwgZGFpbHlfZGVhdGhzKQ0KDQojIENvbnZlcnQgaXQgaW50byBhIHRpbWUgc2VyaWVzDQpkYWlseV9kZWF0aF9oYl96b28gPC0gem9vKHRyZW5kX2RlYXRoX2hiJGRhaWx5X2RlYXRocywgDQogICAgICAgICAgIG9yZGVyLmJ5PWFzLkRhdGUodHJlbmRfZGVhdGhfaGIkZGF0ZSwgZm9ybWF0PSclbS8lZC8lWScpKQ0KDQojIENvbnZlcnQgaXQgaW50byBhIHRpbWUgc2VyaWVzDQpkYWlseV9kZWF0aF9oYl90aW1lc2VyaWVzIDwtICB0aW1lU2VyaWVzOjphcy50aW1lU2VyaWVzKGRhaWx5X2RlYXRoX2hiX3pvbykNCmBgYA0KDQojIyBTdGVwIDEgOiBWaXN1YWxpemUgdGhlIHRpbWUgc2VyaWVzDQoNCmBgYHtyfQ0Kb3JpZ2luYWxfc2VyaWVzX2RlYXRoPC1hdXRvcGxvdChkYWlseV9kZWF0aF9oYl90aW1lc2VyaWVzLCB0cy5jb2xvdXIgPSAnIzVhYjRhYycpKw0KICB4bGFiKCJNb250aCIpICsgDQogIHlsYWIoIlBhdGllbnQgZGllZCIpKw0KICAjc2NhbGVfeF9kYXRlKGJyZWFrcyA9ICIxIG1vbnRoIiwgZGF0ZV9sYWJlbHMgPSAiJWIgLSAleSIgKSsNCiAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA0NSwgdmp1c3QgPSAxLCBoanVzdD0xKSkrDQogIGdndGl0bGUoIlRyZW5kIG9uIERlYXRocyIpICsNCiAgY29sb3JfdGhlbWUoKQ0KDQpnZ3Bsb3RseShvcmlnaW5hbF9zZXJpZXNfZGVhdGgpDQpgYGANCg0KIyMgU3RlcCAyIDogSWRlbnRpZmljYXRpb24gb2YgbW9kZWwgOiAoRmluZGluZyBkOikNCg0KSWRlbnRpZnkgd2hldGhlciB0aGUgdGltZSBzZXJpZXMgaXMgc3RhdGlvbmFyeSAvIG5vbiBzdGF0aW9uYXJ5IHdlIGNhbiB1c2UgQURGIEF1Z21lbnRlZCBEaWNrZXktRnVsbGVyIHRlc3QNCg0KYGBge3J9DQphZGZfdGVzdF9kZWF0aCA8LSBhZGYudGVzdChkYWlseV9kZWF0aF9oYl90aW1lc2VyaWVzKQ0KYWRmX3Rlc3RfZGVhdGgNCmBgYA0KDQpUaGUgdGltZSBzZXJpZXMgaXMgbm90IHN0YXRpb25hcnkgc2luY2Ugd2UgaGF2ZSBhIGhpZ2ggcC12YWx1ZSAocC12YWx1ZSBtdXN0IGJlIFw8IDAuMDUpLg0KU28gd2UgYXBwbHkgZGlmZmVyZW5jZQ0KDQpgYGB7cn0NCmZpcnN0X2RpZmZfZGVhdGg8LSBkaWZmKGRhaWx5X2RlYXRoX2hiX3RpbWVzZXJpZXMpDQphZGZfdGVzdDFfZGVhdGggPC0gYWRmLnRlc3QobmEub21pdChmaXJzdF9kaWZmX2RlYXRoKSkNCmFkZl90ZXN0MV9kZWF0aA0KYGBgDQoNCkNyZWF0ZSBhIGRhdGFmcmFtZSB0byBjb21wYXJlDQoNCmBgYHtyfQ0KYWRmX2RhdGFfZGVhdGggPC0gZGF0YS5mcmFtZShEYXRhID0gYygiT3JpZ2luYWwiLCAiRmlyc3QtT3JkZXJlZCIpLA0KICAgICAgICAgICAgICAgICAgICAgICBEaWNrZXlfRnVsbGVyID0gYyhhZGZfdGVzdF9kZWF0aCRzdGF0aXN0aWMsIGFkZl90ZXN0MV9kZWF0aCRzdGF0aXN0aWMpLA0KICAgICAgICAgICAgICAgICAgICAgICBwX3ZhbHVlID0gYyhhZGZfdGVzdF9kZWF0aCRwLnZhbHVlLGFkZl90ZXN0MV9kZWF0aCRwLnZhbHVlKSkNCmFkZl9kYXRhX2RlYXRoDQpgYGANCg0KSW5pdGlhbGx5IHRoZSBwLXZhbHVlIGlzIGhpZ2ggd2hpY2ggaW5kaWNhdGVzIHRoYXQgdGhlIFRpbWUgU2VyaWVzIGlzIG5vdCBzdGF0aW9uYXJ5Lg0KU28gd2UgYXBwbHkgZGlmZmVyZW5jZSAxIHRpbWUuDQpBZnRlciB0aGUgZmlyc3QgZGlmZmVyZW5jZSwgdGhlIHAtdmFsdWUgXDwgc2lnbmlmaWNhbmNlIGxldmVsICgwLjA1KSBTbyB3ZSBjYW4gY29uY2x1ZGUgdGhhdCB0aGUgZGlmZmVyZW5jZSBkYXRhIGFyZSBzdGF0aW9uYXJ5Lg0KU28gZGlmZmVyZW5jZSAoZCA9IDEpDQoNCk90aGVyIG1ldGhvZDoNCg0KYGBge3J9DQpuZGlmZnMoZGFpbHlfZGVhdGhfaGJfdGltZXNlcmllcykNCmBgYA0KDQpMZXQncyBwbG90IHRoZSBGaXJzdCBPcmRlciBEaWZmZXJlbmNlIFNlcmllcw0KDQpPcmRlciBvZiBmaXJzdCBkaWZmZXJlbmNlDQoNCmBgYHtyfQ0KDQpmaXJzdF9kaWZmX2RlYXRoPC0gZGlmZihkYWlseV9kZWF0aF9oYl90aW1lc2VyaWVzKQ0KcDwtIGF1dG9wbG90KGZpcnN0X2RpZmZfZGVhdGgsIHRzLmNvbG91ciA9ICcjNWFiNGFjJykgKw0KICB4bGFiKCJNb250aCIpICsgDQogIHlsYWIoIkRFQVRIIikrDQogIyBzY2FsZV94X2RhdGUoYnJlYWtzID0gIjEgbW9udGgiLCBkYXRlX2xhYmVscyA9ICIlYiAtICV5IiApKw0KICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDQ1LCB2anVzdCA9IDEsIGhqdXN0PTEpKSsNCiAgZ2d0aXRsZSgiRmlyc3QtT3JkZXIgRGlmZmVyZW5jZSBTZXJpZXMiKSArDQogIGNvbG9yX3RoZW1lKCkNCg0KZ2dwbG90bHkocCkNCmBgYA0KDQojIyBTdGVwIDMgRXN0aW1hdGUgdGhlIHBhcmFtZXRlcnMgKEZpbmRpbmcgcCBhbmQgcSkNCg0KRm9yIG91ciBtb2RlbCBBUklNQSAocCxkLHEpLCB3ZSBmb3VuZCBkID0gMSwgdGhlIG5leHQgc3RlcCBpcyB0byBnZXQgdGhlIHZhbHVlcyBvZiBwIGFuZCBxLCB0aGUgb3JkZXIgb2YgQVIgYW5kIE1BIHBhcnQuDQpQbG90IEFDRiBhbmQgUEFDRiBjaGFydHMgdG8gaWRlbnRpZnkgcSBhbmQgcCByZXNwZWN0aXZlbHkuDQoNCmBgYHtyfQ0KcGFyKG1mcm93PWMoMiwyKSkNCmFjZl9kZWF0aCAgPC0gYWNmMShmaXJzdF9kaWZmX2RlYXRoLCBjb2w9Mjo3LCBsd2Q9NCkNCnBhY2ZfZGVhdGggPC0gYWNmMShmaXJzdF9kaWZmX2RlYXRoLCAgcGFjZiA9IFRSVUUsIGNvbD0yOjcsIGx3ZD00KQ0KYGBgDQoNClRoZSBBQ0YgYW5kIFBBQ0YgcGxvdHMgdGFwZXJzIHRvIHplbyBhdCBvbmUgcG9pbnQuDQpTbyBpdCBmb2xsb3dzIGFuIEFSTUEgbW9kZWwNCg0KVGhlIFBBQ0Ygbm9uLXplcm8gdmFsdWUgYXQgbGFnIDMgU28gdGhlIGRhdGEgbWF5IGZvbGxvdyBhbiBBUigzKSBtb2RlbA0KDQpUaGUgQUNGIG5vbiAtIHplcm8gdmFsdWUgYXQgbGFnIDIgU28gaXQgY2FuIGJlIE1BICgyKQ0KDQpTbyB3ZSBwcm9wb3NlIHRocmVlIEFSTUEgbW9kZWxzIGZvciB0aGUgZGlmZmVyZW5jZWQgZGF0YTogQVJNQShwLHEpIEFSTUEoMywyKSwgQVJNQSgzLDApIGFuZCBBUk1BKDAsMikuDQoNClRoYXQgaXMsIGZvciB0aGUgb3JpZ2luYWwgdGltZSBzZXJpZXMsIHdlIHByb3Bvc2UgdGhyZWUgQVJJTUEgbW9kZWxzLEFSSU1BKHAsZCxxKSBBUklNQSgzLDEsMiksIEFSSU1BKDMsMSwwKSBhbmQgQVJNQSgwLDEsMikuDQoNCiMjIFN0ZXAgNCBCdWlsZCB0aGUgQVJJTUEgbW9kZWwNCg0KIyMjIE1hbnVhbCBBUklNQQ0KDQpgYGB7cn0NCmFyaW1hX2ZpdF9kdGhfMSA9IEFyaW1hKGRhaWx5X2RlYXRoX2hiX3RpbWVzZXJpZXMsIG9yZGVyID0gYygzLDEsMikpDQphcmltYV9maXRfZHRoXzIgPSBBcmltYShkYWlseV9kZWF0aF9oYl90aW1lc2VyaWVzLCBvcmRlciA9IGMoMywxLDApKQ0KYXJpbWFfZml0X2R0aF8zID0gQXJpbWEoZGFpbHlfZGVhdGhfaGJfdGltZXNlcmllcywgb3JkZXIgPSBjKDAsMSwyKSkNCmBgYA0KDQpgYGB7cn0NCnN1bW1hcnkoYXJpbWFfZml0X2R0aF8xKQ0Kc3VtbWFyeShhcmltYV9maXRfZHRoXzIpDQpzdW1tYXJ5KGFyaW1hX2ZpdF9kdGhfMykNCmBgYA0KDQpBbm90aGVyIHdheSBvZiBjaGVja2luZyBBSUMNCg0KYGBge3J9DQp0ZXhyZWc6OnNjcmVlbnJlZyhsaXN0KGFyaW1hX2ZpdF9kdGhfMSwgYXJpbWFfZml0X2R0aF8yLCBhcmltYV9maXRfZHRoXzMpLA0KICAgICAgICAgICAgICAgIGN1c3RvbS5tb2RlbC5uYW1lcyA9YygiQVJJTUEoMywxLDIpIiwiQVJJTUEoMywxLDApIiwiQVJJTUEoMCwxLDIpIiksDQogICAgICAgICAgICAgICAgY2VudGVyID0gVFJVRSwNCiAgICAgICAgICAgICAgICB0YWJsZSA9IEZBTFNFKQ0KYGBgDQoNCkJhc2VkIG9uIHRoaXMsIEFSSU1BIG1vZGVsICgzLDEsMikgc2VlbXMgYmVzdC4NCldlIGNhbiB2ZXJpZnkgdGhlIHNhbWUgdXNpbmcgYXV0b21hdGVkIG1ldGhvZA0KDQojIyMgQXV0b21hdGVkIEFSSU1BDQoNCmBgYHtyfQ0KYXV0b19hcmltYV9maXRfZGVhdGggPC0gYXV0by5hcmltYShsYWcoZGFpbHlfZGVhdGhfaGJfdGltZXNlcmllcyksDQogICAgICAgICAgICAgICAgICBzZWFzb25hbD1GQUxTRSwNCiAgICAgICAgICAgICAgICAgIHN0ZXB3aXNlPUZBTFNFLA0KICAgICAgICAgICAgICAgICAgYXBwcm94aW1hdGlvbj1GQUxTRSwNCiAgICAgICAgICAgICAgICAgIHRyYWNlID0gVFJVRQ0KICAgICAgICAgICAgICAgICAgKQ0KYXV0b19hcmltYV9maXRfZGVhdGgNCmBgYA0KDQpIb3dldmVyIEF1dG9tYXRlZCBBUklNQSBhbHNvIGNvbmZpcm1zIHRoYXQgdGhlIEFSSU1BKDIsMSwzKSBzZWVtcyBnb29kIGJhc2VkIG9uIEFJQw0KDQpgYGB7cn0NCmNvZWZfZHRoPC1sbXRlc3Q6OmNvZWZ0ZXN0KGF1dG9fYXJpbWFfZml0X2RlYXRoKQ0KY29lZl9kdGgNCmBgYA0KDQoqKk1vZGVsIFNlbGVjdGlvbiBDcml0ZXJpYSA6KioNCg0KQVJJTUEgbW9kZWxzIHdpdGggbWluaW11bSBBSUMsIFJNU0UgYW5kIE1BUEUgY3JpdGVyaWEgd2VyZSBjaG9zZW4gYXMgdGhlIGJlc3QgbW9kZWxzLg0KQmFzZWQgb24gQWthaWtlIEluZm9ybWF0aW9uIENyaXRlcmlvbiAoQUlDKSBhYm92ZSwgYW4gQVJJTUEoMiwgMSwgMykgbW9kZWwgc2VlbXMgYmVzdC4NCg0KIyMgU3RlcCA1IENoZWNrIGZvciBEaWFnbm9zdGljcw0KDQpMZXQncyBwbG90IHRoZSBkaWFnbm9zdGljcyB3aXRoIHRoZSByZXN1bHRzIHRvIG1ha2Ugc3VyZSB0aGUgbm9ybWFsaXR5IGFuZCBjb3JyZWxhdGlvbiBhc3N1bXB0aW9ucyBmb3IgdGhlIG1vZGVsIGhvbGQuDQpJZiB0aGUgcmVzaWR1YWxzIGxvb2sgbGlrZSB3aGl0ZSBub2lzZSwgcHJvY2VlZCB3aXRoIGZvcmVjYXN0IGFuZCBwcmVkaWN0aW9uLCBvdGhlcndpc2UgcmVwZWF0IHRoZSBtb2RlbCBidWlsZGluZy4NCg0KYGBge3J9DQpyZXNfZHRoIDwtY2hlY2tyZXNpZHVhbHMoYXV0b19hcmltYV9maXRfZGVhdGgsIHRoZW1lID0gY29sb3JfdGhlbWUoKSkNCnJlc19kdGgNCmBgYA0KDQpUaGUgQUNGIHBsb3Qgb2YgdGhlIHJlc2lkdWFscyBmcm9tIHRoZSBBUklNQSgyLDEsMykgbW9kZWwgc2hvd3MgdGhhdCBhbGwgYXV0byBjb3JyZWxhdGlvbnMgYXJlIHdpdGhpbiB0aGUgdGhyZXNob2xkIGxpbWl0cyBleGNlcHQgMSwgQnV0IHBvcnRtYW50ZWF1IHRlc3QgKExqdW5nLUJveCB0ZXN0KSBzaG93cyBhIGhpZ2hlciBwLXZhbHVlIHdoaWNoIG1lYW5zIHRoZSB2YWx1ZXMgYXJlIGluZHBlbmRlbnQuDQppLmUgV2UgZmFpbCB0byByZWplY3QgdGhlIG51bGwgaHlwb3RoZXNpcw0KDQpGaXR0aW5nIHRoZSBBUklNQSBtb2RlbCB3aXRoIHRoZSBleGlzdGluZyBkYXRhDQoNCkxldCdzIHBsb3QgdGhlIGFjdHVhbHMgYWdhaW5zdCB0aGUgZml0dGVkIHZhbHVlcw0KDQoqKkNvbnZlcnQgbW9kZWwgYW5kIHRpbWUgc2VyaWVzIHRvIGRhdGFmcmFtZSBmb3IgcGxvdHRpbmcqKg0KDQpgYGB7cn0NCmRhaWx5X2RlYXRoX2hiX3RpbWVzZXJpZXNfZGF0YSA8LSBmb3J0aWZ5KGRhaWx5X2RlYXRoX2hiX3RpbWVzZXJpZXMpICU+JSANCiAgY2xlYW5fbmFtZXMoKSAlPiUgDQogIHJlbW92ZV9yb3duYW1lcyAlPiUgDQogIHJlbmFtZSAoZGF0ZSA9IGluZGV4LA0KICAgICAgICAgIGRlYXRoID0gZGF0YSklPiUgDQogIG11dGF0ZShpbmRleCA9IHNlcSgxOm5yb3coZGFpbHlfZGVhdGhfaGJfdGltZXNlcmllcykpKQ0KICANCmFyaW1hX2ZpdF9kdGhfcmVzaWQgPC0gdHMoZGFpbHlfZGVhdGhfaGJfdGltZXNlcmllc1sxOm5yb3coZGFpbHlfZGVhdGhfaGJfdGltZXNlcmllcyldKSAtIHJlc2lkKGF1dG9fYXJpbWFfZml0X2RlYXRoKQ0KDQphcmltYV9maXRfZHRoX2RhdGEgPC0gZm9ydGlmeShhcmltYV9maXRfZHRoX3Jlc2lkKSAlPiUgDQogIGNsZWFuX25hbWVzKCkgJT4lIA0KICBtdXRhdGUoZGF0YSA9IHJvdW5kKGRhdGEsMikpJT4lIA0KICBtdXRhdGUgKGRhdGEgPSBpZmVsc2UoZGF0YSA8IDAsMCxkYXRhKSkNCg0KZml0X2V4aXN0aW5nX2R0aF9kYXRhIDwtIGRhaWx5X2RlYXRoX2hiX3RpbWVzZXJpZXNfZGF0YSAlPiUgDQogIGlubmVyX2pvaW4oYXJpbWFfZml0X2R0aF9kYXRhLCBieSA9IGMoImluZGV4IikpDQoNCg0KYGBgDQoNCioqUGxvdHRpbmcgdGhlIHNlcmllcyBhbG9uZyB3aXRoIHRoZSBmaXR0ZWQgdmFsdWVzKioNCg0KYGBge3J9DQoNCmZpdF9leGlzdGluZ19kdGhfcGxvdCA8LSBmaXRfZXhpc3RpbmdfZHRoX2RhdGEgJT4lIA0KICAgbXV0YXRlIChkYXRhID0gaWZlbHNlKGRhdGEgPCAwLDAsZGF0YSkpICU+JSANCiAgbXV0YXRlKGRhdGUgPSBhcy5EYXRlKGRhdGUpKSAlPiUgDQogIGdncGxvdCgpKw0KICBhZXMoeD1kYXRlLCB5ID0gZGVhdGgpKw0KICBnZW9tX2xpbmUoY29sb3IgPSIjNWFiNGFjIikrDQogIGdlb21fbGluZShhZXMoeD0gZGF0ZSwgeSA9IGRhdGEpLCBjb2xvdXIgPSAicmVkIiApKw0KICB4bGFiKCJNb250aCIpICsgDQogIHlsYWIoIkRlYXRocyByZXBvcnRlZCIpKw0KICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDQ1LCB2anVzdCA9IDEsIGhqdXN0PTEpKSsNCiAgc2NhbGVfeF9kYXRlKGJyZWFrcyA9ICIxIG1vbnRoIiwgZGF0ZV9sYWJlbHMgPSAiJWIgLSAleSIgKSsNCiAgZ2d0aXRsZSgiRml0dGluZyB0aGUgQVJJTUEgbW9kZWwgd2l0aCBleGlzdGluZyBkYXRhIikgKw0KICBjb2xvcl90aGVtZSgpDQoNCmdncGxvdGx5KGZpdF9leGlzdGluZ19kdGhfcGxvdCkNCg0KYGBgDQoNCiMjIFN0ZXAgNiBGb3JlY2FzdCB1c2luZyB0aGUgbW9kZWwNCg0KKipEYXRhIFByZXBhcmF0aW9uIDoqKg0KDQpgYGB7cn0NCmZvcmVjYXN0X2R0aF9tb2RlbCA8LSBmb3JlY2FzdChhdXRvX2FyaW1hX2ZpdF9kZWF0aCxsZXZlbCA9IGMoODAsIDk1KSwgaCA9IDMwKSANCg0KI0NvbnZlcnQgdGhlIG1vZGVsIHRvIGRhdGFmcmFtZSBmb3IgcGxvdHRpbmcNCg0KIyBOZWdhdGl2ZSB2YWx1ZXMgb2YgdGhlIENJIGludGVydmFsIGFyZSBjb25zaWRlcmVkIGFzIDANCg0KZm9yZWNhc3RfZHRoX21vZGVsX2RhdGEgPC0gZm9ydGlmeShmb3JlY2FzdF9kdGhfbW9kZWwpICU+JSANCiAgY2xlYW5fbmFtZXMoKSAlPiUgDQogIG11dGF0ZShkYXRhID0gcm91bmQoZGF0YSwyKSwNCiAgICAgICAgIGZpdHRlZD0gcm91bmQoZml0dGVkLDIpKSAgJT4lIA0KICBtdXRhdGUgKGxvXzgwID0gaWZlbHNlKGxvXzgwIDwgMCwwLGxvXzgwKSwNCiAgICAgICAgICBsb185NSA9IGlmZWxzZShsb185NSA8IDAsMCxsb185NSkNCiAgICAgICAgICApDQoNCmZvcmVjYXN0X3N0YXJ0X2RhdGUgPC0gYXMuRGF0ZShtYXgoZGFpbHlfZGVhdGhfaGJfdGltZXNlcmllc19kYXRhJGRhdGUpKzEpDQpmb3JlY2FzdF9lbmRfZGF0ZSA8LSBhcy5EYXRlKGZvcmVjYXN0X3N0YXJ0X2RhdGUrMjkpDQoNCmZvcmVjYXN0X2R0aF9kYXRhIDwtIGZvcmVjYXN0X2R0aF9tb2RlbF9kYXRhICU+JSANCiAgZmlsdGVyKCEoaXMubmEocG9pbnRfZm9yZWNhc3QpKSkgJT4lIA0KICBtdXRhdGUoZGF0ZSA9IHNlcShmb3JlY2FzdF9zdGFydF9kYXRlLGZvcmVjYXN0X2VuZF9kYXRlLCBieSA9MSkpICU+JSANCnNlbGVjdCgtZGF0YSwtZml0dGVkLCAtaW5kZXgpICANCg0KZml0dGVkX2R0aF9kYXRhIDwtIGZvcmVjYXN0X2R0aF9tb2RlbF9kYXRhICU+JSANCiAgZmlsdGVyKCEoaXMubmEoZGF0YSkpKSAlPiUgDQogIGlubmVyX2pvaW4oZGFpbHlfZGVhdGhfaGJfdGltZXNlcmllc19kYXRhLCBieSA9IGMoImluZGV4IikpICU+JSANCiAgbXV0YXRlKGRhdGUgPSBhcy5EYXRlKGRhdGUpKSAlPiUgDQpzZWxlY3QoZGF0ZSwgZGF0YSwgZml0dGVkKSANCg0KYGBgDQoNCioqUGxvdHRpbmcgdGhlIFZhY2NpbmF0aW9uIHNlcmllcyBwbHVzIHRoZSBmb3JlY2FzdCBhbmQgODAgLSA5NSUgcHJlZGljdGlvbiBpbnRlcnZhbHMqKg0KDQpgYGB7cn0NCg0KYW5ub3RhdGlvbiA8LSBkYXRhLmZyYW1lKA0KICAgeCA9IGMoYXMuRGF0ZSgiMTMtMDYtMjAyMSIsIiVkLSVtLSVZIiksYXMuRGF0ZSgiMTEtMTAtMjAyMSIsIiVkLSVtLSVZIikpLA0KICAgeSA9IGMoMzAsNDApLA0KICAgbGFiZWwgPSBjKCJQQVNUIiwgIkZVVFVSRSIpDQopDQoNCiNUaW1lIHNlcmllcyBwbG90cyBmb3IgdGhlIG5leHQgNjAgZGF5cyBhY2NvcmRpbmcgdG8gYmVzdCBBUklNQSBtb2RlbHMgd2l0aCA4MCXigJM5NSUgQ0kuDQpmb3JlY2FzdF9kYXRhX2R0aF9wbG90IDwtIGZpdHRlZF9kdGhfZGF0YSAlPiUgDQogIGdncGxvdCgpKw0KICBnZW9tX2xpbmUoYWVzKHg9IGRhdGUsIHkgPSBkYXRhKSwgY29sb3IgPSAiIzVhYjRhYyIpKw0KICAjZ2VvbV9saW5lKGFlcyh4PSBkYXRlLCB5ID0gZml0dGVkKSwgY29sb3VyID0gInJlZCIgKSsNCiAgZ2VvbV9saW5lKGFlcyh4PSBkYXRlLCB5ID1wb2ludF9mb3JlY2FzdCksIGNvbG9yID0iYmx1ZSIsIGRhdGEgPSBmb3JlY2FzdF9kdGhfZGF0YSApKw0KICBnZW9tX3JpYmJvbihhZXMoeCA9IGRhdGUsIHkgPSBwb2ludF9mb3JlY2FzdCwgeW1pbiA9IGxvXzgwLCB5bWF4ID0gaGlfODApLCANCiAgICAgICAgICAgICAgZGF0YSA9IGZvcmVjYXN0X2R0aF9kYXRhLCBhbHBoYSA9IDAuMywgZmlsbCA9ICJncmVlbiIpKw0KICBnZW9tX3JpYmJvbihhZXMoeCA9IGRhdGUsIHkgPSBwb2ludF9mb3JlY2FzdCwgeW1pbiA9IGxvXzk1LCB5bWF4ID0gaGlfOTUpLCANCiAgICAgICAgICAgICAgZGF0YSA9IGZvcmVjYXN0X2R0aF9kYXRhLCBhbHBoYSA9IDAuMSkrDQogIGdlb21fdmxpbmUoYWVzKHhpbnRlcmNlcHQ9YXMubnVtZXJpYyhtYXgoZGF0ZSkpKSxjb2xvcj0iI2YxYTM0MCIsIGxpbmV0eXBlPSJkYXNoZWQiLGRhdGEgPSBmaXR0ZWRfZHRoX2RhdGEpKw0KICBnZ3RpdGxlKCJQcm9qZWN0aW9uIG9mIG5ldyBEZWF0aHMiKSArDQogIHhsYWIoIk1vbnRoIikgKyANCiAgeWxhYigiRGVhdGggcmVwb3J0ZWQiKSsNCiAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA0NSwgdmp1c3QgPSAxLCBoanVzdD0xKSkrDQogIGNvbG9yX3RoZW1lKCkrDQogIHNjYWxlX3hfZGF0ZShicmVha3MgPSAiMSBtb250aCIsIGRhdGVfbGFiZWxzID0gIiViIC0gJXkiICkrDQogIGdlb21fdGV4dChkYXRhPWFubm90YXRpb24sIA0KICAgICAgICAgICAgYWVzKCB4PXgsIHk9eSwgbGFiZWw9bGFiZWwpLCAgICAgICAgICAgICAgICAgIA0KICAgICAgICAgICAgY29sb3I9ImJsdWUiLCANCiAgICAgICAgICAgIHNpemU9NCApDQogICANCg0KZ2dwbG90bHkoZm9yZWNhc3RfZGF0YV9kdGhfcGxvdCApDQpgYGANCg==